Stability of hexagonal solidification patterns 
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We investigate the dynamics of cellular solidification patterns using three-dimensional phase-field 
simulations. The cells can organize into stable hexagonal patterns or exhibit unsteady evolutions. 
We identify the relevant secondary instabilities of regular hexagonal arrays and find that the stability 
boundaries depend significantly on the strength of crystalline anisotropy. We also find multiplet 
states that can be reached by applying well-defined perturbations to a pre-existing hexagonal array. 



I. INTRODUCTION 



Hexagonal patterns arise in many non-equilibrium systems when a two-dimensional translational symmetry is broken 
by a dynamic instability 0. Well-known examples are Benard-Marangoni convection or Faraday surface waves. In 
directional solidification, where alloy samples are pulled with a fixed speed V in an externally imposed temperature 
gradient G, the planar solidification front becomes unstable above a critical value of the ratio V/G. Close above this 
threshold, the front develops shallow cells that often form regular hexagonal arrays nil; as V/G increases, cells 
deepen and eventually form dendrites with lateral branches. Besides the interest of this sequence of morphological 
transitions as a model system for the physics of pattern formation 0|, it also has practical importance because the 
morphology of the solidification front determines the spatial distribution of alloy components and defects in the 
solidified material that, in turn, influence its mechanical properties. 

We focus here on two issues: the influence of crystalline anisotropy on pattern selection, and the existence of 
multiplet states consisting of several symmetry-broken cells. These questions have been examined in two dimensions 
(2D), both in thin-sample experiments |5|,^ JJ and numerically by boundary integral and phase-field methods 0,0,1^. 
The situation is less clear in three dimensions (3D), where the analysis of experimental data must account for the 
presence of convection in the bulk liquid, and numerical investigations have so far been restricted to rather short 
length and time scales Hol [Tp . Here, we use an efficient formulation of the phase-field method [Ti| and a multi- 
scale simulation algorithm [l3j to perform extensive quantitative 3D simulations of a symmetric model of directional 
solidification. 

The influence of crystallographic effects on solidification patterns was already noticed in early experimental studies 
1^1^. It is by now well established that for freely growing dendrites, pattern selection is governed by small anisotropics 
of the solid-liquid surface tension and interface kinetics. More surprisingly, it was recently shown that anisotropy plays 
a crucial role even close above the instability threshold in directional solidification, where the cell shapes are round 
For thin-sample experiments, where the system is quasi-two-dimensional, the anisotropy of the one-dimensional 
interface is controlled by the orientation of the crystal with respect to the sample plane ^6] . Stable and regular steady 
states are obtained only for sufficiently strong anisotropy; for weak anisotropy, cells may split in two or disappear, 
and the solidification front never reaches steady state. In numerical simulations without anisotropy, cells are stable 
only in a narrow parameter range close to the onset of the primary instability; when anisotropy is included, stable 
solutions exist even far from threshold . Here, we investigate in detail the stability of three-dimensional hexagonal 
cells and identify the relevant secondary instabilities. We find that the anisotropy has the same strong effect on array 
stability as in 2D. 

Furthermore, we find triplet solutions; that is, groups of three cells that grow with their tips close together. This 
is the 3D equivalent of doublet fingers found in 2D both in experiments and simulations |M S 13 ■ Similar 3D growth 
structures have been found in simulations of free growth in a channel |^ , and as transients in experiments . In 
our simulations, triplets can be generated in a controlled manner starting from a pre-existing hexagonal array. 



II. MODEL AND IMPLEMENTATION 



An alloy of overall composition cq is pulled with velocity y in a linear temperature field with its gradient G 
directed along the z axis. In the frame of the sample, the temperature is given by T = Tq + G{z — Vt), where 
solidification starts at t = from a flat equilibrium interface of temperature Tq at z = 0. This is the so-called frozen 
temperature approximation, valid for slow solidiflcation and similar thermal conductivities of both phases. The phase 
field (j) distinguishes between solid {4> = 1) and liquid (0 = —1). The concentration jump between solid and liquid 
Ac is assumed constant. Instead of the concentration field, which exhibits rapid variations across the interface, the 
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TABLE L Simulation parameters. For all simulations, Wo = 1, D = 2, the grid spacing is Aa::=0.8 and the time step At = 0.05. 
To obtain a "physical" capillary anisotropy 64, a slightly higher value £4 has to be used in the model to correct for grid effects 
(see Ref. O). 

£4 €4 a To 5 

0.00 0.0029 3.482 1.018 -0.0180 

0.01 0.0129 3.407 0.9965 0.00345 

0.02 0.0231 3.330 0.975 0.0255 

0.03 0.0334 3.262 0.955 0.0472 



dimcnsionless field u ~ (c— co)/Ac— (1 + (j))/2 is used. It is constant for an equilibrium interface and hence equivalent 
to a chemical potential. The time evolution of the fields is given by 

T{n)dt(f>^ — — , (1) 

dtu = DV'^u + 9(0/2, (2) 

where fi = V(/)/|V(/)| is the unit normal to the interface, T{fi) is an orientation-dependent relaxation time, D is the 
solute diffusivity, and 5F / 5<j) denotes the functional derivative of the Lyapounov functional 

F - 1 + fw + (. - ^) , (3) 

where W{n) is the orientation-dependent interface thickness. It — mAc/G is the thermal length with m being the 
liquidus slope in the phase diagram, and /((/>) — — 0^/2 + 0''/4 is a double- well potential. The quantity u— {z — Vt) /It 
is the local supersaturation, and its product with g(0) — (p— 20'^/3 -I- (p^/b introduces a free energy difference between 
the two phases by tilting the double well, which provides the driving force for the phase transformation. Anisotropic 
surface tension with cubic symmetry and one principal axis aligned with the growth direction is implemented by 
W{n) = Woa{n) with a{h) = {1 - 3€4)[{l + 4:e4{nl + n'^ + nj)/ {1 - 364)]- In the thin-interface limit [3, this model is 
equivalent to the classic free-boundary problem of solidification, in which u satisfies the diffusion equation in solid and 
liquid, the Stefan condition at the interface, Vn = n- [Z?Vii|s — _DVu|;], where Vn is the normal velocity of the interface, 
and the gradients of u are evaluated on the solid and liquid side of the interface, and the generalized Gibbs-Thomson 
condition at the interface, 

Mint = --^ doy^ [a{n) + dg.a{h)] - (3{h)Vn, (4) 

It , 
i—i 

where 9i are the local angles between h and the two local principal directions on the interface, Ki are the principal 
curvatures, do = aiWo/a is the capillary length, and P{n) — aiT/{Wa)[l — a2aW^ / {Dt)] is the kinetic coefficient, 
with oi = 5a/2/8 and 02 = 0.6267. Taking into account grid corrections as explained in detail in Ref. the interface 
kinetics can be eliminated {f3 = 0) by setting T{n) — to(1 — 3S)[{l + 4:6{n^ + ny + n'^)/{l — 36)] and by choosing specific 
values for the constants tq, (5, and a that depend on D, the capillary anisotropy £4 and the discretization; the values 
used for our simulations are listed in tabled 

The involved physical length scales are the thermal length It, the diffusion length Id = D/V, and the capillary 
length do. The stability of a planar interface is controlled by the parameter v = It /Id- In experiments, v is usually 
varied by changing V at fixed G; in the simulations, it is more efficient to vary G, and hence It at fixed do/lo- 
Typical experimental values for dt^/ljj are of the order 10^^. However, due to computational constraints, such small 
values are too costly to simulate. Indeed, convergence of the phase-field method requires W/Id <C 1 and Vt/W < 1. 
The latter, combined with the requirement = (that is, Dt/W''^ ^ a ^ W/do), yields W/Id < \/do/lD- Since 
the discretization is on the scale of W, whereas typical cell sizes are of order Id or larger, the required number of 
grid points increases when do/'n decreases. For the value do /Id = 2.5 x 10~^ used in most of our simulations, the 
morphological instability occurs at Vc = 3.75 (instead of i^c = 2 in the limit do/lo 0). 

The above model has parallel liquidus and solidus lines and equal diffusivities in the solid and the liquid, whereas 
alloys exhibit strongly different diffusivities in the two phases and a temperature-dependent concentration jump. 
Nevertheless, this model has been shown to reproduce even complex qualitative features of pattern selection that are 
seen in 2D experiments 3i ^-^d computationally remains the most efficient model that captures the salient features 
of directional solidification, despite recent advances on more realistic alloy models |15| . 
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FIG. 1: Evolution of cellular patterns for v = 12.5 and (a) £4 = 0.03 and (b) £4 = 0.02. Left: top views of the growth front at 
tD/l^ = 320, 640, and 1996.8 (the greyscale is proportional to the surface height; box size: Lx/Id = Ly/lo ~ 12.8); right: 3D 
view of the interface at the end of the run. 

III. INFLUENCE OF ANISOTROPY 

Simulations begin with an unstable planar steady-state interface, and numerical noise triggers the morphological 
instability. The subsequent evolution of the front dramatically depends on the anisotropy strength. For an anisotropy 
of £4 = 0.03, some of the initial cells are eliminated, but most survive, become stable, and slowly drift over small 
distances until they form a perfectly regular lattice (Fig. QJi). In other runs, hexagonal arrays with defects are 
formed; in all cases, the interface shape converges to a stable steady-state solution. In contrast, for a slightly lower 
anisotropy (£4 = 0.02) with otherwise identical conditions, the cells do not reach steady-state, and cell-elimination 
and cell-splitting events occur up to the end of our simulation runs (Fig. ^). 

As a first step toward the understanding of pattern selection, we investigate the stability of regular hexagonal arrays. 
We start from initial conditions as described in Fig. [51 and allow the cells to reach steady state at moderate values of 
v. Then, we increase v in small steps until a secondary instability threshold is reached. Two different instabilities are 
encountered. For large spacings, the cells oscillate and form a superlattice pattern shown in Fig. |21 Three equivalent 
sublattices of hexagonal symmetry emerge, with -y/S times larger spacing than the basic pattern. All cells on each 
sublattice oscillate in phase with an exponential growth in amplitude over time. Eventually, the cells tip-split, and 
the subsequent dynamics is unsteady and similar to Fig. ^p. For smaller spacings, cell elimination occurs: all cells 
located on the same sublattice fall behind the others and are overgrown. The subsequent evolution is again unsteady. 
Finally, for v very close to the onset of the primary instability, a few examples of isolated cell eliminations occurred 
that are probably the signature of an Eckhaus instability. 

The geometry of the instability modes is entirely determined by the underlying hexagonal symmetry of the initial 
pattern . Similar superlattice oscillations have been reported in simulations of two-dimensional amplitude equations 
|l7t IT^ , and experimentally observed in fluid systems [l^, '2^ . The new feature here is the presence of anisotropy. 
Figure shows the stability boundaries of oscillatory and cell elimination instabilities computed for several values of 
£4. Whereas the cell elimination mode is only weakly sensitive to anisotropy, the critical cell spacing for the onset of 
the oscillatory mode is strongly anisotropy-dependent. Without anisotropy, the stability boundaries of the oscillatory 
and cell elimination modes cross and no stable states were found for ly > 9. For £4 — 0.01, a narrow range of stable 
states was found to extend to the largest values of ly investigated; this stable band becomes larger with increasing 
anisotropy. These findings correlate well with the overall behavior of the large-scale runs of Fig. depicted in the 
inset of Fig. |21 Stable cellular arrays appear spontaneously only when a sufhciently large band of stable states is 
available. 

Remarkably, the structure of the stability diagram is identical to the one obtained for 2D in Ref. (Qj. The cell 
elimination and cell oscillation modes found here for hexagonal arrays are the equivalent of the well-known steady and 
oscillatory period-doubling instabilities found in 2D systems. To further investigate this similarity, we have performed 
two-dimensional simulations of our phase-field model and found that, for fixed simulation parameters, the stability 
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FIG. 2: Left: sketch of the simulated geometry: in a box of lateral size x Ly with periodic boundary conditions, m rows of n 
cells may be generated by starting from a flat surface perturbed by the superposition of three plane waves; we define a — L^jn 
and 6 = Ly/m. Right: snapshot of an oscillatory instability that develops spontaneously in a box with 24 cells (m = 4, n = 6) 
for v = 5.4, a/lo = 3.307, b/lo = 2.880, and £4 = 0.01. 
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FIG. 3: Main plot: stability boundaries for hexagonal arrays for different anisotropies: 0% (circles), 1% (squares) 2% (dia- 
monds), and 3% (triangles). Hexagonal arrays are stable between the limits of cell elimination (dotted lines and open symbols) 
and oscillations (dasli-dottcd linos and full symbols). The full lino is the neutral stability line of the primary instability. Inset: 
behavior of large runs started from a flat surface. Filled squares: stable hexagons, open circles: unsteady evolution, dashed 
line: primary instability threshold. 



diagrams in 2D and 3D arc not identical; however, the critical spacings for the onset of instability differ by less than 
20% in all cases investigated. 



IV. OTHER INSTABILITY MODES 



To investigate how the spatial structure of the oscillatory modes depends on the symmetry of the underlying 
steady-state solution, simulations starting from elongated hexagons were performed. By varying the aspect ratio 
of the simulation box, patterns ranging from rows of elongated hexagons {b/a > \/3/2) over symmetric hexagons 
{b/a = V3/2) to squares (b/a = 1/2) can be generated. For a square lattice, there are only two equivalent sublattices, 
and a two-sublattice oscillation mode emerges. For rows of elongated hexagons, close to a one-dimensional symmetry. 
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FIG. 4: Instabilities of deformed hexagons for v = 5 and £4 = 0. The lengths a and b axe defined in Fig. |21 the dashed line 
corresponds to regular hexagons, the dotted line to squares. 




FIG. 5: Triplets at = 12.5, a/lo ~ 5.333, b/lo = 4.640, and £4 = 0.02. Six cells were simulated and the final picture was 
replicated four times to yield a clearer view. 

in-phase oscillations of entire rows occur. For strongly elongated hexagons, a tilt instability is found that corresponds 
to a breaking of the parity symmetry: the cells become asymmetric and drift along their longer axis. However, these 
states are unstable with respect to subsequent oscillations. Figure 0] shows the symmetry of the fastest instability 
in the plane of the two lengths a and b defined in Fig. [2 We find a balloon of stable arrays consisting of regular or 
slightly elongated hexagons. 



V. TRIPLETS 



To search for multiplet states, we use a procedure inspired by the experimental method used in Ref. 7] to create 
doublets in thin-sample directional solidification. A UV-absorbing laser dye was used as solute, and by illuminating 
the sample with UV light through a mask, the temperature field was perturbed for a limited time with a well-defined 
spatial pattern. We implement a similar procedure by adding to the temperature field a perturbation AT{x,y) that 
has the symmetry of one of the three sublattices of Fig. |21 and whose minima are located over the "holes" between 
three cells. This promotes the growth of the surrounding cells, and their tips approach each other to form triplet 
fingers. Once the perturbation is switched off, the system relaxes to triplet states that remain stable for the duration 
of our simulations for large enough v and large wavelengths (Fig. . For lower z/, the triplets exhibit instabilities 
similar to the usual cells, but no detailed investigation was carried out. Triplets never appeared spontaneously in our 
simulations: only unsteady evolutions resulted from random initial conditions, even for control parameters for which 
stable triplet states exist. 
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VI. CONCLUSION 



We have studied the stabihty of hexagonal arrays and triplet fingers, and found that the crystalline anisotropy plays 
a crucial role in determining the stability bounds. Remarkably, despite the cubic symmetry of the anisotropy, the 
selected patterns are hexagonal, while square arrays are always unstable. Square cellular states have been calculated 
long ago using a numerical scheme valid only for shallow cells pll |. However, the geometry of the used simulation cell 
suppresses all large-scale oscillatory modes, such that no conclusion about their stability can be drawn. Nevertheless, 
the existence of other than hexagonal patterns cannot be ruled out in general, since it is known that the weakly 
nonlinear behavior of solidification fronts strongly depends on the alloy characteristics Q , and we have not carried 
out an exhaustive survey of the (large) parameter space. 

The main approximations made in the present work are the constant concentration jump, the equal solute diffu- 
sivities, and the large value of do /Id- The two first can be easily levied by the use of a more general model [THj l: 
the third remains a serious computational challenge. We expect most of our results to remain qualitatively valid for 
more realistic systems, since they depend only on the the underlying symmetries and the stability of the cell tips. In 
that sense, it is interesting to note that unsteady evolutions similar to Fig. QJd and transient isolated multiplcts have 
recently been observed in 3D experiments on transparent alloys |l4j . 
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